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A granular gas subjected to a permanent injection of energy is described by means of hydrody- 
namic equations derived from a moment expansion method. The method uses as reference function 
not a Maxwellian distribution /m but a distribution /o = $ /m , such that $ adds a fourth cumulant 
k to the velocity distribution. The formalism is applied to a stationary conductive case showing 
that the theory fits extraordinarily well the results coming from our molecular dynamic simulations 
once we determine k as a function of the inelasticity of the particle-particle collisions. The shape of 
k is independent of the size N of the system. 



I. INTRODUCTION 

Granular systems subjected to a sufficiently strong excitation may have a fluid-like behavior plp f. From the very 
beginning several authors have attempted to derive hydrodynamic equations for these systems [|3|]4|7^ If the excitation 
of the granular system is through a permanent injection of energy, the fluid system may stabilize to a low density 
stationary gaseous state which necessarily is a noncquilibrium state and it usually is inhomogeneous as well. To develop 
the basic features of the theory of gaseous granular systems we restrict the analysis to the simplifying inelastic hard 
sphere model (IHS) §. 

Many authors studying granular gases have put particular attention to studying the spontaneous homogeneous 
cooling of a granular system using periodic boundary conditions . This time dependent state is called homogeneous 
cooling state (HCS) and the understanding of its properties has been improving through many articles [f7|-|lo|. A crucial 
breakthrough was the realization by Goldshtein and Shapiro |Q that the homogeneous cooling distribution function 
has a scaling property with respect to the instantaneous temperature. Such distribution — which we will be calling 
/hcs — is known in approximate forms [ljjjisj ]. It is known, among other things, that its fourth cumulant k does not 
vanish and that it has a long velocity tail. 

A nonequilibrium inhomogeneous gaseous system, on the other hand, is described by a distorted distribution function 
typically obtained from Boltzmann's equation expanding the distribution either in gradients of the hydrodynamic 
fields (Chapmann-Enskog's method) jl3| or making a moment expansion (Grad's method) [Q. For normal gases the 
expansion is made about the equilibrium Maxwell's distribution. 

In this article we will assume that a low density nonequilibrium granular system has a local distribution function 
which can be obtained expanding about a distribution /o resembling a /hcs in the sense that it has a significantly 
nonvanishing fourth cumulant. We introduce a reference function (see Eq. (||) below) which is a Maxwellian distorted 
by a factor which incorporates the next non trivial cumulant, a fourth cumulant, to the distribution function and then 
we undertake a perturbative moment expansion a la Grad about /o to solve Boltzmann's equation. Some authors 
have done some calculations in this direction but using the gradient expansion (Chapman-Enskog) method ]l5p . The 
point is that, without the notion of equilibrium, we expect that the reference state, /o in our case, should resemble 
more the homogeneous cooling state than the simple Maxwellian. 

We study a two-dimensional system of hard disks, and the moment expansion — in dimension two — is an 8 moment 
expansion: the number density n(r, i), the velocity field w(r, t), the granular temperature field T(r, t), the pressure 
tensor Pij(f,t) and the heat flux vector field Q(r,t). The dynamic variables are not the components of the pressure 
tensor IP itself but the components of the symmetric traceless part pij where Py = p 5ij + and p is the hydrostatic 
pressure. 
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As it can be seen in Grad's article (TJ] or in |l6| the method yields hydrodynamic equations for all the fields 
mentioned above. In particular, the dynamic equations for pij and Q take the place of what would normally be the 
constitutive (transport) equations of standard hydrodynamics. This last point means that we are not assuming any 
constitutive equations whatsoever, their present counterparts are dynamic equations. 

It is well established that, in the case of the IHS model for granular systems, Boltzmann's equation is modified in 
that the restitution coefficient r enters solely in the gain term of the collision integral and it does so in two forms. 
First the gain term has an overall factor r~ 2 and second, the distribution functions appearing in the gain term depend 
upon the precollision velocities, and these velocities depend on r 

When the modified Boltzmann's equation is used the stemming hydrodynamic equations get factors that depend 
on the inelasticity coefficient q, 

1 - r 

(q — in the perfectly elastic case) except that the mass continuity equation and the momentum balance equation 
remain unchanged since mass and momentum continue being microscopically conserved. 

In the context of Boltzmann's equation a dissipative gas satisfies the ideal gas equation of state 

p = nT (2) 

where the granular temperature T is defined in energy units as the average kinetic energy per particle. If we were 
to consider the Boltzmann-Enskog equation, then the inelasticity coefficient would enter through the Enskog collision 
factor Xt an d the equation of state of a normal gas and a dissipative gas would differ, but in the present context 
Eq. (§) holds. 

Usual moment expansion methods (as Grad's is) are appropriate to describe bulk properties. Wall effects are 
not well described unless higher order momenta are included in the expansion, which are not trivial to handle [ p"7| |. 
Already in normal gases hydrodynamic fields may have discontinuities at walls. In a previous article we gave a kinetic 
description of a one-dimensional granular system with theoretical tools such that our description was correct and 
precise up to the walls pl| but we have not generalized yet that type of formalism to higher dimensions, hence, in 
the present article, we use moment expansions. 

In consequence the formalism we are going to present suffers too of the weakness of moment expansions: it is 
unreliable precisely at the points where the boundary conditions should be imposed forcing us to trade the boundary 
conditions for conditions imposed far from the walls based on the actual behavior of the system according to our 
molecular dynamic simulations. 

Our moment expansion method, explained in detail in Sec. II, uses as reference distribution function a distribution 
/o which differs from a Maxwellian distribution in that it has a nonvanishing fourth cumulant, k. 

The method leads to hydrodynamic equations for low density granular systems that depend parametrically on the 
inelasticity coefficient q and k. We apply this hydrodynamic equations to a stationary and purely conductive case 
as in p9| . The value of the fourth cumulant k, or better, the dependence of the fourth cumulant on q is determined 
directly from molecular dynamic simulations, and it turns out to be independent of the size N of the system. The 
predictions that follow from our formalism agree very well with all our simulational data. 

In Sec. [n] we briefly present the moment expansion method, in Sec. H the hydrodynamic equations that follow are 
given and specialized to a purely conductive case and finally in Sec. |lVptheory and simulational results are compared. 
Final comments arc in Sec. 0. 



II. THE MOMENT EXPANSION METHOD FOR GRANULAR SYSTEMS 



Moment expansion methods can summarily be described as follows. Take a velocity distribution function fo(f, c, t) 
which is considered to be the reference function about which an expansion is going to be made. For normal gases the 
natural choice for /o is a local Maxwellian distribution /m written in terms of the peculiar velocity C = c — v(r,t), 
where v(r,t) is the hydrodynamic velocity. Next a set of orthonormal polynomials on C, H a (C), are built in the 
sense that Ho = 1 and 

J H a (C)H b (C)f (r,C,t)dC = S ab (3) 

The polynomials H a are obtained simply building a base of orthonormal polynomials starting from Hq = 1 and from 
first degree upwards. When / is a Maxwellian the H a are Hermite polynomials but in general they are not. 
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Then a solution of the form 



ffi C,t)=(l + J2 H a (C) R a (r, t) f (r, C, t) 



(4) 



is replaced in Boltzmann's equation. The R a are the moments of / with respect to the H a and they can directly be 
related to the moments associated to c (the velocity field v), to pCiCj (the kinematic pressure tensor Pij) and to 
\p C 2 C (the heat flux vector Q). One could go on but we have used polynomials only containing C,, Cj Cj and C 2 Ci 
as Grad did. 

The following step is to derive integrability conditions multiplying the kinetic equation consecutively by the H n and 
then integrating the equation over C. The idea is to do this up to a given order and drop all contributions coming 
from polynomials of degree higher than a chosen value (up to order 3 in our case) . This gives a set of hydrodynamic 
equations for the different moments. 

A key point is the choice of the reference function /q. The two dimensional Maxwellian /m = 
n 2^ry exp[— to C 2 /(2T)] is privileged as the solution describing the equilibrium state of a normal gas. Since in granular 
systems there is no such thing as equilibrium a next best choice, seems a distorted Maxwellian distribution |ll| 

C 4N 



/o = 

where the dimensionless peculiar velocity is 



f 



-|l-C 2 



c 



TO 



(5) 



(6) 



and T is the granular temperature. The coefficient re is the fourth cumulant of / and it depends on the inelasticity 
coefficient q while the coefficients in front of C 2 and C 4 are derived from requiring that fo is normalized and that 
<fC 2 ) 



fo 



T. The fourth cumulant re in dimension two is 



(C 4 )-2(C 2 f 



(7) 



In the case of the homogeneous cooling state, recent articles have justified explicit forms |ll|Jl2[ | for re in the context 
of distribution functions like fo- Their results are well approximated by 



with 



bi = -2, 



bi+b 2 q 
1 + &3 q ' 



62 ~ 13.619, 



4.5969. 



(8) 



(9) 



The rational expression given in (||) is valid within 2.9% for q < 0.08, r = 1 — 2q = 0.84. Notice that n(q = 0) = 
allowing to recover the elastic case. 

We describe quite satisfactorily nonhomogeneous, nonequilibrium stationary states with a K as h ra), but with 
numerical factors different from those in Eq. (||). The latter values are valid only for the homogeneous cooling state 
while our system is kept in a stationary inhomogeneous state with appropriate boundary conditions. 

For any re the resulting distribution obtained with the method summarized above is 



Px 



P(2 + 
Q x (C 2 



-(C 2 X -C 2 ) + : 

i - 2k)C x 



Pxy 



p(2 + n 



2Q (2 + 5re - re 2 ) 2Q (2 + 5re 



fo 



(10) 



where, p yy — p xx and Qq — y ^ p. The distribution shares with /o the first scalar moments: density, temperature 

and fourth cumulant, for any value of q. If re is chosen to be zero then, in Eq. (f[o|), /o — > /m and becomes the 
usual Grad's distribution. Hence the whole method would be the original method devised by Grad and, if re is chosen 
to be Eq. (||) with coefficient values as in Eq. (||), then f would be what we are calling /hcs- 

Given our ignorance regarding granular gases one could, in principle, accept /m or /hcs as legitimate reference 
functions to make the moment expansion. In the following sections we compare the three formalisms (reference 
functions /m, /hcs an d /o, the latter with a re adjusted to the results) concluding that only the one based on /o gives 
acceptable results for a sufficiently large range of q. In fact they are very good. 
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III. THE HYDRODYNAMIC EQUATIONS 



As already mentioned, the inelasticity coefficient q enters the kinetic equation in two different forms, ft appears 
as a factor = nz^p m the gain term and it appears in the expression for the precollision velocities which are 

part of the argument of the distribution functions appearing in the gain term. When the solution is inserted in 
Boltzmann's equation, q enters in a still third form, precisely through the n coefficient given by Eqs. (^|) and (|Tc|). 
Expanding the collisonal term of Boltzmann's equation in powers of q and k, the moment method yields the following 
hydrodynamics equations for a granular gas 

Dn , s 

— + nV-v = 0, (11) 

Dv 

mn— -nF + V ■ P = 0, (12) 

n^ + V.Q + P:V*=- 2 -^l, (13 ) 
Dt T 

dpij d 1 / dQ l dQj dQ k \ 

^— (vkPij) + r I -q— + -a~r ~ Siinz- + 



dt dxk 2 \ dxj dxi dx k 

dvi | dvj dv s ( dvi dvj f dv r \ B(q) 



Prj dx- r +P "dx~ r - S - Prs dx- r +P {dx- + dx~- 5ij dx- r ) = " — ^ ' (W) 



-5T - + (<Wfc) + 7: ^— Qr + X ^ + X ^— Qfe 

oi ax r 2 ox r 2 ox% 2 cte r 



T 9pfcr | ipkr OT p kr dP rs | 3 p dT _ C(q) ^ 
m dx r m dx r m n dx s m dx k 2 r ' 



where D/ Dt = d/dt + v ■ V , 



1 mT . 
2ap V 7T 

is a characteristic relaxation time, a is the diameter of the particles, and the coefficients A, B and C for the generic 
distribution are 

A(q) = [q(l -q)+ 0(q 5 )}(l + ^K+ , 

B{q) = i + h-h 2 + fefl + §1+ A + + °(" 3 ) , (17) 



c(«) 



(i + |«) 

1 + f g - f g 2 + ^(206 + 1267g) - f§|U 2 + Q(q 3 ) + Q( K 3 ) 
(! + §«-£) 



They depend on g and k, but k depends on q. At least for small q the coefficients A, B and C are positive (g varies 
between and |). The coefficient A in Eq. ( |13|) determines the energy dissipation in the system and consequently it 
vanishes in the elastic limit while B and C tend to 1. In the previous equations both p±j and Pij appear depending 
on which of them gives a more compact expression. 



It is our hope that the above hydrodynamics is valid in a wide variety of situations compatible with gaseous states 
and no clustering but in this article we restrict our study to a hydrostatic case. 

We are going to consider the purely conductive regime, with no external force, F = 0. The system of N disks is 
in a rectangular box of dimension L x x L y , with thermal walls, at y = — ^ L y and at y = \ L y , both at temperature 
T , and periodic boundary conditions in the X direction. The case with no external force is quite simple because the 
system is symmetric with respect to y — > — y and the pressure is uniform. If the system were conservative, as a normal 
gas, this would be a homogeneous system at thermal equilibrium, but since the system is dissipative the temperature 
depends on the coordinate y (T has a minimum at the symmetry axis) and the problem is much less trivial. 
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Since the system is purely conductive there is no velocity field, P xy = 0, Q x = and the other fields depend only 
on the coordinate y. The set of equations (|lll-[l5|) becomes 



_ B ~ A p B + A 

(18) 

where the prime indicates derivative with respect to y. Notice that because there is inelasticity the pressure tensor is 
anisotropic in the sense that P xx ^ P yy . In fact P yy — P xx ~ A ~ q and they do not depend on the coordinate y. 

Next we are going to compare the implications of these equations with our molecular dynamics results. To this 
end we should, in principle, solve these hydrodynamic equations using the boundary conditions associated to the 
particular simulations that we have studied. This is not a straightforward task because, as we have mentioned at the 
end of the introduction, the moment expansion method (behind the previous hydrodynamic equations) does not give 
a good description near boundaries. For example, if we impose that a wall behaves as a stochastic wall at temperature 
To, the observed field T is not expected to take that value near the wall. Later on it will be shown how we tackle this 
problem. 

From Eq. (|l^) it is direct to derive that the temperature field satisfies the equation 

TT" + -T' 2 = k 2 a 2 p 2 , (19) 

e^ A \ ABC . (20) 

3A + 2B v ' 

Since the pressure is uniform, and because of Eq. (0), this equation can also be used as an equation for the inverse of 
the density. 

Before proceeding to solve the equations we adimensionalize the problem defining a coordinate £ = y/L yi (— i < 
£ < i) and we also define 



where 



N 



L X Ly 



n{y) = nn*(^) , 



* = (TT) ' =*Toy/% Q;(0, (21) 

T(y) =T T*(C), p =nT p*. 

Once Eq. jis| ) is solved and converted to a solution for the dimensionless number density n*(£) the result is 



AK K| = * J **-;"*® + J- arctanh >-:^ , (22) 

n I?/ y n max n uiax y n max 

where is the value taken by in the middle of the system, ^ = 0. In what follows it is shown that nj^ ax is 

a simple function of the parameter K defined above. The effect of the boundaries is traded in favor of the observed 
values of n* nax . 

The integral condition J ndxdy = N can be cast as 

1/2 i 

n*dt=-, (23) 
but since d£ = (d^/dn*) dn* the integral condition yields, 

= 1 - tanh 2 K , (24) 

77 * 

max 

where n£ is the value that ri*(£) takes at the two boundaries, n£ = n*(±^). 
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Equation ( p2| ) evaluated at £ = \ yields a different condition over n* h and it follows that 

, _ 1 , sinh(2*Q 



(25) 



max 2 AK 

This is a strong result, it says that the dimensionless number density at the center, n^ ax , is determined by K 2 alone, 

where = N Pa/oi, a = L x /L y is the aspect ratio of the box, and pa — (ir/A)(Na 2 / (L x L y )) is the fraction of area 
occupied by the disks. Since A vanishes in the elastic limit and for small inelasticity coefficient q, A sa q while B»l 
and C w 1 then 

K 2 na 2 qn -qNp A . (27) 
a 

In the quasielastic limit then, the control parameter is basically qN pa- This result, for fixed area density, resembles 
what we obtained in the one dimensional case Jl8| , namely, that the relevant control parameter of the one dimensional 
equations is qN. 

The temperature is T*(£) = p* /n*(£) but since the pressure is uniform one may be tempted to use p* — T b * with 
the value for n* b already derived from the theory and the value T b * imposed in the simulation. This would give a bad 
fit however because, as we have been emphasizing, the formalism is not reliable near the boundaries. Therefore we 
choose for p* the value p* — T* in n^ ax . We know n^ ax from Eq. (p5j), and we take the value for T^ in directly from 
our simulations. It may be said that Tj^ in is a parameter to adjust our results. 

The dimensionless heat flux becomes, 




r y^ = -Ybc~ ~ ^ — ■ (28) 



IV. SIMULATION-THEORY COMPARISON 



In this section we compare the simulational results with the values given by three formalisms which use as reference 
function: (i) /m, (m) /hcs & n d (Hi) fo- We are calling /hcs the function like /o but with k defined with the values 
given in Eq. (^|). 

Since the formalisms differ by terms of higher order in q their predictions are quite similar unless qN is large enough. 
Typically the Maxwellian theory and the one based on /hcs are valid until about qN — 20 and are reasonable until 
qN = 40 while the theory based on fo, with an adjusted k, is valid for values of qN up to 200, and reasonably good 
until about qN w 300. 



A. Simulational setup 



We have performed simulations of a two dimensional system of N — 2300, N = 3600, N = 10000, and N = 19600 
inelastic hard disks inside a L x L box with lateral periodic boundary conditions while the upper and lower walls are 
kept at granular temperature To = 1. The area fraction covered by the disks was chosen to be pa = 0.01 (in which 
case the nonideal corrections to the equation of state are less than 2%) while the qN dissipation parameter ranges 
from qN = 10 up to qN = 400. In the N — 2300 case, the smallest simulated system, the ratio between the mean 
free path and the linear size of the system (Knudsen number) is 0.065 and it is smaller in the other cases. This value 
guarantees that not too close to the walls the fluid has a hydrodynamic behavior. The wall temperature To is imposed 
sorting the velocity of the bouncing particles as if they were coming from a heat bath at T — Tq. 

In every simulation the system was relaxed from an initial condition for a sufficiently long time. After the relaxation 
we measure local time averages of the main moments of the distribution (i.e., n, v, T, pij, Q) inside each one of a 
set of square cells. Taking advantage of the translation invariance in the X direction, it is natural to take horizontal 
averages getting, in this way, smooth vertical profiles for the observed hydrodynamic fields. 



6 



B. Theory and simulation comparison 



In order to compare theory and simulation some analysis is needed because size effects are quite noticeable unless 
the number N of particles is above about 3600. We use expression (E5j) for n* ax as our point of contact between 
theory and simulation and proceed to adjust the coefficients bk in Eq. (§f) so that K takes the values that make theory 
give the observed values for n* max . 

With this aim we first write K as a rational expression K = ao yjq (1 + a\ q)/(l + 0,2 q) and find the values of the at 
so that Eq. (|2^) reproduces the observed values of n^ ax . We have to adjust a in spite of its definition, given under 
Eq. (p6|), because the effective values for the number of particles, the global density and the aspect ratio get distorted 
since there is a layer near the thermalizing walls which does not behave hydrodynamically. Once this expression for K 
is fixed, we invert ( p6j ) to obtain values for k as a function of q. The size effect is in a alone and n(q) is approximately 
the same for systems with N larger than about 3600, as shown in Fig. [I]. In this figure there is a solid line which is 
Eq.(|) with 

61 = -22.541, 6 2 = 6.19187, 63 = 79.0362 . (29) 

The discrepancies between this curve and the empirical values of k are less than 2% for the whole range of q considered. 
More in detail, Fig. |l| shows the behavior of n(q) for systems of different size. It is seen that in the case of N = 2300 
(solid circles) n(q) is dependent on the system's size, while the predictions in the cases N = 3600, 10000, and 19600 
differ among themselves by less than 2%. Only the smallest simulated system (N = 2300) departs from this otherwise 
universal shape. 

In the following we present the results corresponding to the N = 10000 case. 

a. The density: As a first step and in order to check the validity of the kinetic description (no clustering, gj) 
we have plotted the final configurational positions of the particles (not shown here) for different values of qN up to 
qN = 400, founding that clustering begins at about qN w 300. In fact, we have detected that the number of collisions 
per unit time increases dramatically shortly after qN — 300. 

As a second step we check the validity of Eq. ( p5| ) (which for fixed geometry depends only on qN) with simulations. 
The top Fig. || shows both n^ ax (growing curves with qN) and n* h (the decreasing curves). It is seen how well the 
values of n* mm (solid circles) come from the fo distribution (there should be no surprise as these are the fitted data), 
compared with the prediction using the other two distributions. This graph shows that n^ ax with the other two 
formalisms is good only up to qN sa 40. 

Also at top Fig. |^ shows the values predicted by Eq. ( p4] ) and the observed values of n* b (see caption), and it is 
seen that all theoretical schemes give values close to the simulational results in the considered range of qN. This is 
the only case where, for large qN, predictions coming from /m are better that those coming from fo, but we do not 
believe there is anything deep here since our moment expansion method is not reliable near walls. 

Figure |^, bottom, compares theory and simulational density profiles for qN = 30 and qN = 200. As it has already 
been explained, the value n^ ax in Eq. ( |22] ) is fixed by Eq. ( Poj ) and there are no extra parameters to adjust. It is 
seen that theory in all cases (/m, /hcs and fo) give good agreement for low values of qN. However predictions for 
qN = 200, when using /m or /hcs fail. It is amazing how large qN can be when f is used. In the last case the 
parameter K is adjusted only from the knowledge of n^ ax and it accurately predicts the behavior in almost all the 
volume. 

From the figure it can be appreciated that near the walls (£ = ±0.5), as mentioned in the introduction, the theory 
does not predict well the behavior of «*(£)• 

b. The temperature: It has already been mentioned that the temperature profiles exhibit a minimum at the center 
and there is a temperature jump at the boundaries. In Fig. || we show the simulational results for T^ in = T*(£ = 0) 
and temperature T b * = T*(£ = ±1/2). It is seen that for qN = 40 the temperature of the fluid by the walls is about 
40% lower than the imposed value. This effect is due to dissipation and in ID it has been shown to be a Q(qN) 
effect @. 

Because at present we have no theory to describe the temperature jump that takes place near the thermal walls, 
and we know that Grad's method does not give good results near boundaries, we have chosen the observed value of 
the temperature at mid height, Tj^ in , as the value to use in the formalism. The observed values decrease with qN and, 
in the case N = 10000, we have adjusted them with the following expression 

ln(T min ) = -0.0163259 qN + 5.27656 10" 5 {qN) 2 

-1.86396 10~ 7 (qN) 3 + 4.241 10~ 10 (qN) 4 

-4.23878 10" 13 (qN) 5 (30) 

whose faithfulness is shown in Fig. 
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Since the pressure is uniform the temperature profile can be written directly as 



T*. 77* 

rTi*//-\ mm max \ 

Figure |] shows the simulational and theoretical temperature profiles for some values of qN. The three upper curves 
show the profiles for small qN values (qN — 10, 20, 30) and it is seen that the predictions using fo (solid line) 
give an excellent fit to the simulational data while for the results obtained using /m or /hcs (dashed and light lines 
respectively) the fit is only fair. The lowest curves shows T profiles for larger values of qN (qN = 50, 100, 200). In this 
case only the formalism based on / give an acceptable description and the agreement is very good up to qN — 200. 
In the case of qN = 275, not shown, there is still a reasonable agreement when fo is used. 

c. The pressure: Now that the values of n* max and T^ in are known and since the pressure is uniform (both in theory 
and simulationally) and theory asserts that p* — T^ m n max , then we have the value of p* . 

From the value for the pressure and the set of equations ( |l8|) one gets the theoretical values for P yy in terms of the 
pressure tensor. These values are compared with our simulational data in Fig. ^| Again the comparison is good when 
the formalism with fo is used and it is only fairly good with the other two formalisms. 

As it has already been mentioned after Eq. (|l8|), the two diagonal terms of the pressure tensor are not equal because 
the system is anisotropic. 

d. The heat current: Equation ( pSj ) gives the theoretical expression for the heat current Q* . At top Fig. ^| shows 
simulational data and theoretical predictions for small values of qN. The three formalism predict well the observed 
values. 

At bottom in Fig. || it is possible to see that for qN = 200 only the formalism using fo fits the observed data, while 
the other two fail badly. In the case of qN = 275 the agreement is still reasonable within about 3%. 

V. FINAL COMMENTS 

In this article we have studied a bidimensional granular system of N inelastic hard disks (normal restitution 
coefficient r). The system is placed in a rectangular box and it is kept in a stationary regime with upper and lower 
walls at granular temperature To. The lateral walls are periodic. The quantity qN has been used as control parameter, 
where q = and simulations were made with different values of N ranging from N — 2300 to 19600. If the system 
were conservative it would remain in a perfectly homogeneous state at temperature To with a homogeneous area 
density pa = 0.01. Because there is dissipation the dimensionless number density n* has a maximum in the middle, 
n max , and the dimensionless temperature T* has a minimum, T^ in , also at the center of the system. The pressure is 
uniform and p* = < lax T* in . 

Hydrodynamic equations were derived using a moment expansion method. This method has the fourth cumulant, 
re, of the velocity distribution as a parameter and three possible re's where considered. These are: re = which implies 
that the reference distribution function in a Maxwellian, and two re's of the form (|§|), one with parameters bk as in 
(^) which corresponds to using /hcs as the reference function and finally using the values bk, Eq. (]2^), numerically 
determined to ensure that the correct values of the density at the middle of the system come out. The last case gives 
our reference function fo. 

The empirical rational form of re = re(g), Eq. (||), in the case of the successful distribution fo turns out to be 
independent of the size of the system. 

When this fo is used, the comparison between the predictions and simulation results for n max , p 

*, p xx, Py V against 

dissipation, and the density, temperature and heat flux profiles are very good. 

The obvious conclusion is that the formalism based on fo gives an excellent description of the behavior of the 
system for values of qN up to 200 and it gives a reasonable description up to qN nearly 300 (slightly above qN ~ 300 
clustering begins), while the other formalisms fail beyond about qN = 40. 

Even though the fourth cumulant re of fo is obtained to fit the results of a particular hydrodynamic regime, we 
expect that the theory with this expression for re is still valid for other regimes. 
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FIG. 1. The fourth cumulant k against q for systems of different size. The empirical values of ft are represented by triangles 
for N = 19600, empty circles for N = 10000, empty squares for N = 3600, and solid circles for N = 2300. The solid line 
corresponds to our empirical fit. See text. 



FIG. 2. At top the predicted and observed values for the density n^ ax at the center of the box (£ = 0) and nl near 
the boundaries (£ = ±0.5). The solid (open) circles correspond to the simulational values for n* ax (n£), the light-dashed 
(heavy-dashed) line corresponds to the theoretical prediction using /hcs (/m). The solid line corresponds to our empirical 
adjustment (see text). At bottom the predicted and observed density profiles for two different values of qN. The open (solid) 
circles correspond to qN — 30 (qN = 200). The light-dashed (heavy-dashed) line corresponds to the theoretical prediction 
using /hcs (/m). The solid line corresponds to the prediction stemming from fo- 

FIG. 3. Temperature T* in at the center of the channel (open circles) against dissipation qN. The solid line corresponds to a 
fit using Eq. (|30J) . The solid circles are the observed values for the temperature T 6 * near the boundaries 

FIG. 4. From top to bottom temperature profiles corresponding to qN = 5, 10, 30, 50, 100 and 200. The empty circles 
are the simulational results, the dashed line (light solid line) correspond to the predictions obtained with /m (/hcs)- The heavy 
solid lines are the theoretical results when fo is used. 

FIG. 5. Simulational and theoretical values for the components of the pressure tensor against dissipation qN. At top (bottom) 
the P xx (P yy ) component. In light-solid line (dashed line) the theoretical results when using /hcs distribution (/m distribution). 
The heavy-solid line is the theoretical result when using the fo distribution. 

FIG. 6. Simulational and theoretical values for the heat flux profiles. At top the three curves and set of simulational data 
correspond to qN — 5 (circles), qN = 20 (squares) and qN = 30 (rhombus). The dashed and light-solid lines correspond to 
the predictions using /m and /hcs respectively, the heavy solid line corresponds to fo- At bottom are the observed values of Q y 
when qN = 200 and, the theoretical predictions when using fo (heavy solid line), /m (dashed line) and /hcs (light-solid line). 
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